{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Set-up instructions:** this notebook give a tutorial on the forecasting learning task supported by `sktime`.\n",
    "On binder, this should run out-of-the-box.\n",
    "\n",
    "To run this notebook as intended, ensure that `sktime` with basic dependency requirements is installed in your python environment.\n",
    "\n",
    "To run this notebook with a local development version of sktime, either uncomment and run the below, or `pip install -e` a local clone of the `sktime` `main` branch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from os import sys\n",
    "# sys.path.append(\"..\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forecasting with sktime\n",
    "\n",
    "In forecasting, past data is used to make temporal forward predictions of a time series. This is notably different from tabular prediction tasks supported by `scikit-learn` and similar libraries.\n",
    "\n",
    "\n",
    "<img src=\"img/forecasting.png\" width=750 />\n",
    "\n",
    "`sktime` provides a common, `scikit-learn`-like interface to a variety of classical and ML-style forecasting algorithms, together with tools for building pipelines and composite machine learning models, including temporal tuning schemes, or reductions such as walk-forward application of `scikit-learn` regressors.\n",
    "\n",
    "**Section 1** provides an overview of common forecasting workflows supported by `sktime`.\n",
    "\n",
    "**Section 2** discusses the families of forecasters available in `sktime`.\n",
    "\n",
    "**Section 3** discusses advanced composition patterns, including pipeline building, reduction, tuning, ensembling, and autoML.\n",
    "\n",
    "**Section 4** gives an introduction to how to write custom estimators compliant with the `sktime` interface.\n",
    "\n",
    "Further references:\n",
    "* for further details on how forecasting is different from supervised prediction à la `scikit-learn`, and pitfalls of misdiagnosing forecasting as supervised prediction, have a look at [this notebook](./01a_forecasting_sklearn.ipynb)\n",
    "* for a scientific reference, take a look at [our paper on forecasting with sktime](https://arxiv.org/abs/2005.08067) in which we discuss `sktime`'s forecasting module in more detail and use it to replicate and extend the M4 study."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "\n",
    "* [1. Basic forecasting workflows](#chapter1)\n",
    "    * [1.1 Data container format](#section_1_1)\n",
    "    * [1.2 Basic deployment workflow - batch fitting and forecasting](#section_1_2)\n",
    "        * [1.2.1 Basic deployment workflow in a nutshell](#section_1_2_1)\n",
    "        * [1.2.2 Forecasters that require the horizon already in `fit`](#section_1_2_2)\n",
    "        * [1.2.3 Forecasters that can make use of exogeneous data](#section_1_2_3)\n",
    "        * [1.2.4 Prediction intervals](#section_1_2_4)      \n",
    "    * [1.3 basic evaluation workflow - evaluating a batch of forecasts against ground truth observations](#section_1_3)   \n",
    "        * [1.3.1 The basic batch forecast evaluation workflow in a nutshell - function metric interface](#section_1_3_1)\n",
    "        * [1.3.2 The basic batch forecast evaluation workflow in a nutshell - metric class interface](#section_1_3_2)           \n",
    "    * [1.4 advanced deployment workflow: rolling updates & forecasts](#section_1_4) \n",
    "        * [1.4.1 updating a forecaster with the update method](#section_1_4_1)    \n",
    "        * [1.4.2 moving the \"now\" state without updating the model](#section_1_4_2)   \n",
    "        * [1.4.3 walk-forward predictions on a batch of data](#section_1_4_3)  \n",
    "    * [1.5 advanced evaluation worfklow: rolling re-sampling and aggregate errors, rolling back-testing](#section_1_5)         \n",
    "* [2. Forecasters in sktime - main families](#chapter2)\n",
    "    * [2.1 exponential smoothing, theta forecaster, autoETS from statsmodels](#section_2_1)\n",
    "    * [2.2 ARIMA and autoARIMA](#section_2_2)\n",
    "    * [2.3 BATS and TBATS](#section_2_3)    \n",
    "    * [2.4 Facebook prophet](#section_2_4)  \n",
    "    * [2.5 State Space Model (Structural Time Series)](#section_2_5)  \n",
    "* [3. Advanced composition patterns - pipelines, reduction, autoML, and more](#chapter3)\n",
    "    * [3.1 Reduction: from forecasting to regression](#section_3_1)\n",
    "    * [3.2 Pipelining, detrending and deseasonalization](#section_3_2)    \n",
    "        * [3.2.1 The basic forecasting pipeline](#section_3_2_1)\n",
    "        * [3.2.2 The Detrender as pipeline component](#section_3_2_2) \n",
    "        * [3.2.3 Complex pipeline composites and parameter inspection](#section_3_2_3)        \n",
    "    * [3.3 Parameter tuning](#section_3_3)      \n",
    "        * [3.3.1 Basic tuning using ForecastingGridSearchCV](#section_3_3_1)  \n",
    "        * [3.3.2 Tuning of complex composites](#section_3_3_2)        \n",
    "        * [3.3.3 Selecting the metric and retrieving scores](#section_3_3_3) \n",
    "    * [3.4 autoML aka automated model selection, ensembling and hedging](#section_3_4) \n",
    "        * [3.4.1 autoML aka automatic model selection, using tuning plus multiplexer](#section_3_4_1)   \n",
    "        * [3.4.2 autoML: selecting transformer combinations via OptimalPassthrough](#section_3_4_2)  \n",
    "        * [3.4.3 simple ensembling strategies](#section_3_4_3)   \n",
    "        * [3.4.4 Prediction weighted ensembles and hedge ensembles](#section_3_4_4)  \n",
    "* [4. Extension guide - implementing your own forecaster](#chapter4)        \n",
    "* [5. Summary](#chapter5)          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### package imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Basic forecasting workflows <a class=\"anchor\" id=\"chapter1\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section explains the basic forecasting workflows, and key interface points for it.\n",
    "\n",
    "We cover the following three workflows:\n",
    "\n",
    "* basic deployment workflow: batch fitting and forecasting\n",
    "* basic evaluation workflow: evaluating a batch of forecasts against ground truth observations\n",
    "* advanced deployment workflow: fitting and rolling updates/forecasts\n",
    "* advanced evaluation worfklow: using rolling forecast splits and computing split-wise and aggregate errors, including common back-testing schemes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Data contanier format<a class=\"anchor\" id=\"section_1_1\"></a>\n",
    "\n",
    "All workflows make common assumptions on the input data format.\n",
    "\n",
    "`sktime` uses `pandas` for representing time series:\n",
    "\n",
    "* `pd.Series` for univariate time series and sequences\n",
    "* `pd.DataFrame` for multivariate time series and sequences\n",
    "\n",
    "The `Series.index` and `DataFrame.index` are used for representing the time series or sequence index. `sktime` supports pandas integer, period and timestamp indices.\n",
    "\n",
    "NOTE: at current time (v0.6x), forecasting of multivariate time seres is not a stable functionality, this is a priority roadmap item. Multivariate exogeneous time series are part of stable functionality."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example:** as the running example in this tutorial, we use a textbook data set, the Box-Jenkins airline data set, which consists of the number of monthly totals of international airline passengers, from 1949 - 1960. Values are in thousands. See \"Makridakis, Wheelwright and Hyndman (1998) Forecasting: methods and applications\", exercises sections 2 and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.datasets import load_airline\n",
    "from sktime.utils.plotting import plot_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:05.849641Z",
     "iopub.status.busy": "2021-04-10T16:07:05.849188Z",
     "iopub.status.idle": "2021-04-10T16:07:06.049023Z",
     "shell.execute_reply": "2021-04-10T16:07:06.049499Z"
    }
   },
   "outputs": [],
   "source": [
    "y = load_airline()\n",
    "\n",
    "# plotting for visualization\n",
    "plot_series(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, users are expected to use the in-built loading functionality of `pandas` and `pandas`-compatible packages to load data sets for forecasting, such as `read_csv` or the `Series` or `DataFrame` constructors if data is available in another in-memory format, e.g., `numpy.array`.\n",
    "\n",
    "`sktime` forecasters may accept input in `pandas`-adjacent formats, but will produce outputs in, and attempt to coerce inputs to, `pandas` formats.\n",
    "\n",
    "NOTE: if your favourite format is not properly converted or coerced, kindly consider to contribute that functionality to `sktime`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Basic deployment workflow - batch fitting and forecasting<a class=\"anchor\" id=\"section_1_2\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest use case workflow is batch fitting and forecasting, i.e., fitting a forecasting model to one batch of past data, then asking for forecasts at time point in the future.\n",
    "\n",
    "The steps in this workflow are as follows:\n",
    "\n",
    "1. preparation of the data\n",
    "2. specification of the time points for which forecasts are requested. This uses a `numpy.array` or the `ForecastingHorizon` object.\n",
    "3. specification and instantiation of the forecaster. This follows a `scikit-learn`-like syntax; forecaster objects follow the familiar `scikit-learn` `BaseEstimator` interface.\n",
    "4. fitting the forecaster to the data, using the forecaster's `fit` method\n",
    "5. making a forecast, using the forecaster's `predict` method\n",
    "\n",
    "The below first outlines the vanilla variant of the basic deployment workflow, step-by-step.\n",
    "\n",
    "At the end, one-cell workflows are provided, with common deviations from the pattern (Sections 1.2.1 and following)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 1 - preparation of the data\n",
    "\n",
    "as discussed in Section 1.1, the data is assumed to be in `pd.Series` or `pd.DataFrame` format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.datasets import load_airline\n",
    "from sktime.utils.plotting import plot_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in the example, we use the airline data set.\n",
    "y = load_airline()\n",
    "plot_series(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 2 - specifying the forecasting horizon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to specify the forecasting horizon and pass that to our forecasting algorithm.\n",
    "\n",
    "There are two main ways:\n",
    "\n",
    "* using a `numpy.array` of integers. This assumes either integer index or periodic index (`PeriodIndex`) in the time series; the integer indicates the number of time points or periods ahead we want to make a forecast for. E.g., `1` means forecast the next period, `2` the second next period, and so on.\n",
    "* using a `ForecastingHorizon` object. This can be used to define forecast horizons, using any supported index type as an argument. No periodic index is assumed.\n",
    "\n",
    "Forecasting horizons can be absolute, i.e., referencing specific time points in the future, or relative, i.e., referencing time differences to the present. As a default, the present is that latest time point seen in any `y` passed to the forecaster.\n",
    "\n",
    "`numpy.array` based forecasting horizons are always relative; `ForecastingHorizon` objects can be both relative and absolute. In particular, absolute forecasting horizons can only be specified using `ForecastingHorizon`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### using a numpy forecasting horizon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:06.222860Z",
     "iopub.status.busy": "2021-04-10T16:07:06.222322Z",
     "iopub.status.idle": "2021-04-10T16:07:06.224682Z",
     "shell.execute_reply": "2021-04-10T16:07:06.225165Z"
    }
   },
   "outputs": [],
   "source": [
    "fh = np.arange(1, 37)\n",
    "fh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will ask for monthly predictions for the next three years, since the original series period is 1 month.\n",
    "In another example, to predict only the second and fifth month ahead, one could write:\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "fh = np.array([2, 5])  # 2nd and 5th step ahead\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Using a `ForecastingHorizon` based forecasting horizon\n",
    "\n",
    "The `ForecastingHorizon` object takes absolute indices as input, but considers the input absolute or relative depending on the `is_relative` flag.\n",
    "\n",
    "`ForecastingHorizon` will automatically assume a relative horizon if temporal difference types from `pandas` are passed; if value types from `pandas` are passed, it will assume an absolute horizon.\n",
    "\n",
    "To define an absolute `ForecastingHorizon` in our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.base import ForecastingHorizon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:06.229051Z",
     "iopub.status.busy": "2021-04-10T16:07:06.228401Z",
     "iopub.status.idle": "2021-04-10T16:07:06.230685Z",
     "shell.execute_reply": "2021-04-10T16:07:06.231173Z"
    }
   },
   "outputs": [],
   "source": [
    "fh = ForecastingHorizon(\n",
    "    pd.PeriodIndex(pd.date_range(\"1961-01\", periods=36, freq=\"M\")), is_relative=False\n",
    ")\n",
    "fh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`ForecastingHorizon`-s can be converted from relative to absolute and back via the `to_relative` and `to_absolute` methods. Both of these conversions require a compatible `cutoff` to be passed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoff = pd.Period(\"1960-12\", freq=\"M\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fh.to_relative(cutoff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fh.to_absolute(cutoff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 3 - specifying the forecasting algorithm\n",
    "\n",
    "To make forecasts, a forecasting algorithm needs to be specified. This is done using a `scikit-learn`-like interface. Most importantly, all `sktime` forecasters follow the same interface, so the preceding and remaining steps are the same, no matter which forecaster is being chosen.\n",
    "\n",
    "For this example, we choose the naive forecasting method of predicting the last seen value. More complex specifications are possible, using pipeline and reduction construction syntax; this will be covered later in Section 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.naive import NaiveForecaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = NaiveForecaster(strategy=\"last\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 4 - fitting the forecaster to the seen data\n",
    "\n",
    "Now the forecaster needs to be fitted to the seen data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster.fit(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 5 - requesting forecasts\n",
    "\n",
    "Finally, we request forecasts for the specified forecasting horizon. This needs to be done after fitting the forecaster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = forecaster.predict(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting predictions and past data\n",
    "plot_series(y, y_pred, labels=[\"y\", \"y_pred\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2.1 the basic deployment workflow in a nutshell<a class=\"anchor\" id=\"section_1_2_1\"></a>\n",
    "\n",
    "for convenience, we present the basic deployment workflow in one cell.\n",
    "This uses the same data, but different forecaster: predicting the latest value observed in the same month."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.base import ForecastingHorizon\n",
    "from sktime.forecasting.naive import NaiveForecaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 1: data specification\n",
    "y = load_airline()\n",
    "\n",
    "# step 2: specifying forecasting horizon\n",
    "fh = np.arange(1, 37)\n",
    "\n",
    "# step 3: specifying the forecasting algorithm\n",
    "forecaster = NaiveForecaster(strategy=\"last\", sp=12)\n",
    "\n",
    "# step 4: fitting the forecaster\n",
    "forecaster.fit(y)\n",
    "\n",
    "# step 5: querying predictions\n",
    "y_pred = forecaster.predict(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional: plotting predictions and past data\n",
    "plot_series(y, y_pred, labels=[\"y\", \"y_pred\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2.2 forecasters that require the horizon already in `fit` <a class=\"anchor\" id=\"section_1_2_2\"></a>\n",
    "\n",
    "Some forecasters need the forecasting horizon provided already in `fit`. Such forecasters will produce informative error messages when it is not passed in `fit`. All forecaster will remember the horizon when already passed in `fit` for prediction. The modified workflow to allow for such forecasters in addition is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 1: data specification\n",
    "y = load_airline()\n",
    "\n",
    "# step 2: specifying forecasting horizon\n",
    "fh = np.arange(1, 37)\n",
    "\n",
    "# step 3: specifying the forecasting algorithm\n",
    "forecaster = NaiveForecaster(strategy=\"last\", sp=12)\n",
    "\n",
    "# step 4: fitting the forecaster\n",
    "forecaster.fit(y, fh=fh)\n",
    "\n",
    "# step 5: querying predictions\n",
    "y_pred = forecaster.predict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:06.475031Z",
     "iopub.status.busy": "2021-04-10T16:07:06.473831Z",
     "iopub.status.idle": "2021-04-10T16:07:06.613175Z",
     "shell.execute_reply": "2021-04-10T16:07:06.613700Z"
    }
   },
   "source": [
    "#### 1.2.3 forecasters that can make use of exogeneous data<a class=\"anchor\" id=\"section_1_2_3\"></a>\n",
    "\n",
    "Many forecasters can make use of exogeneous time series, i.e., other time series that are not forecast, but are useful for forecasting `y`. Exogeneous time series are always passed as an `X` argument, in `fit`, `predict`, and other methods (see below). Exogeneous time series should always be passed as `pandas.DataFrames`. Most forecasters that can deal with exogeneous time series will assume that the time indices of `X` passed to `fit` are a super-set of the time indices in `y` passed to `fit`; and that the time indices of `X` passed to `predict` are a super-set of time indices in `fh`, although this is not a general interface restriction. Forecasters that do not make use of exogeneous time series still accept the argument (and do not use it internally).\n",
    "\n",
    "The general workflow for passing exogeneous data is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 1: data specification\n",
    "y = load_airline()\n",
    "# we create some dummy exogeneous data\n",
    "X = pd.DataFrame(index=y.index)\n",
    "\n",
    "# step 2: specifying forecasting horizon\n",
    "fh = np.arange(1, 37)\n",
    "\n",
    "# step 3: specifying the forecasting algorithm\n",
    "forecaster = NaiveForecaster(strategy=\"last\", sp=12)\n",
    "\n",
    "# step 4: fitting the forecaster\n",
    "forecaster.fit(y, X=X, fh=fh)\n",
    "\n",
    "# step 5: querying predictions\n",
    "y_pred = forecaster.predict(X=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: as in workflows 1.2.1 and 1.2.2, some forecasters that use exogeneous variables may also require the forecasting horizon only in `predict`. Such forecasters may also be called with steps 4 and 5 being\n",
    "```python\n",
    "forecaster.fit(y, X=X)\n",
    "y_pred = forecaster.predict(fh=fh, X=X)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2.4 prediction intervals<a class=\"anchor\" id=\"section_1_2_4\"></a>\n",
    "\n",
    "`sktime` provides a unified interface to return prediction interval when forecasting. This is possible directly in the `predict` function, by setting the `return_pred_int` argument to `True`. The `predict` method then returns a second argument, Not all forecasters are capable of returning prediction intervals, in which case an error will be raised.\n",
    "\n",
    "Obtaining prediction intervals can be done as part of any workflow involving `predict`, by adding the argument `return_pred_int` - below, we illustrate this by modifying the basic workflow in Section 1.2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.theta import ThetaForecaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple workflow\n",
    "y = load_airline()\n",
    "fh = np.arange(1, 13)\n",
    "\n",
    "forecaster = ThetaForecaster(sp=12)\n",
    "forecaster.fit(y)\n",
    "\n",
    "# setting return_pred_int argument to True; alpha determines percentiles\n",
    "#  intervals are lower = alpha/2-percentile, upper = (1-alpha/2)-percentile\n",
    "alpha = 0.05  # 2.5%/97.5% prediction intervals\n",
    "y_pred, y_pred_ints = forecaster.predict(fh, return_pred_int=True, alpha=alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`y_pred_ints` is a `pandas.DataFrame` with columns `lower` and `upper`, and rows the indices for which forecasts were made (same as in `y_pred`). Entries are lower/upper (as column name) bound of the nominal alpha predictive interval for the index in the same row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_ints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pretty-plotting the predictive interval forecasts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plot_series(y, y_pred, labels=[\"y\", \"y_pred\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plot_series(y, y_pred, labels=[\"y\", \"y_pred\"])\n",
    "ax.fill_between(\n",
    "    ax.get_lines()[-1].get_xdata(),\n",
    "    y_pred_ints[\"lower\"],\n",
    "    y_pred_ints[\"upper\"],\n",
    "    alpha=0.2,\n",
    "    color=ax.get_lines()[-1].get_c(),\n",
    "    label=f\"{1 - alpha}% prediction intervals\",\n",
    ")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: this should be turned into a one-liner, by moving this to `utils.plotting` - contributions are appreciated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 basic evaluation workflow - evaluating a batch of forecasts against ground truth observations<a class=\"anchor\" id=\"section_1_3\"></a>\n",
    "\n",
    "It is good practice to evaluate statistical performance of a forecaster before deploying it, and regularly re-evaluate performance if in continuous deployment. The evaluation workflow for the basic batch forecasting task, as solved by the workflow in Section 1.2, consists of comparing batch forecasts with actuals. This is sometimes called (batch-wise) backtesting.\n",
    "\n",
    "The basic evaluation workflow is as follows:\n",
    "\n",
    "1. splitting a representatively chosen historical series into a temporal training and test set. The test set should be temporally in the future of the training set.\n",
    "2. obtaining batch forecasts, as in Section 1.2, by fitting a forecaster to the training set, and querying predictions for the test set\n",
    "3. specifying a quantitative performance metric to compare the actual test set against predictions\n",
    "4. computing the quantitative performance on the test set\n",
    "5. testing whether this performance is statistically better than a chosen baseline performance\n",
    "\n",
    "NOTE: step 5 (testing) is currently not supported in `sktime`, but is on the development roadmap. For the time being, it is advised to use custom implementations of appropriate methods (e.g., Diebold-Mariano test; stationary confidence intervals).\n",
    "\n",
    "NOTE: note that this evaluation set-up determines how well a given algorithm would have performed on past data. Results are only insofar representative as future performance can be assumed to mirror past performance. This can be argued under certain assumptions (e.g., stationarity), but will in general be false. Monitoring of forecasting performance is hence advised in case an algorithm is applied multiple times."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example:** In the example, we will us the same airline data as in Section 1.2. But, instead of predicting the next 3 years, we hold out the last 3 years of the airline data (below: `y_test`), and see how the forecaster would have performed three years ago, when asked to forecast the most recent 3 years (below: `y_pred`), from the years before (below: `y_train`). \"how\" is measured by a quantitative performance metric (below: `mean_absolute_percentage_error`). This is then considered as an indication of how well the forecaster would perform in the coming 3 years (what was done in Section 1.2). This may or may not be a stretch depending on statistical assumptions and data properties (caution: it often is a stretch - past performance is in general not indicative of future performance)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 1 - splitting a historical data set in to a temporal train and test batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.model_selection import temporal_train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = load_airline()\n",
    "y_train, y_test = temporal_train_test_split(y, test_size=36)\n",
    "# we will try to forecast y_test from y_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting for illustration\n",
    "plot_series(y_train, y_test, labels=[\"y_train\", \"y_test\"])\n",
    "print(y_train.shape[0], y_test.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 2 - making forecasts for y_test from y_train\n",
    "\n",
    "This is almost verbatim the workflow in Section 1.2, using `y_train` to predict the indices of `y_test`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can simply take the indices from `y_test` where they already are stored\n",
    "fh = ForecastingHorizon(y_test.index, is_relative=False)\n",
    "\n",
    "forecaster = NaiveForecaster(strategy=\"last\", sp=12)\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "\n",
    "# y_pred will contain the predictions\n",
    "y_pred = forecaster.predict(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting for illustration\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### steps 3 and 4 - specifying a forecasting metric, evaluating on the test set\n",
    "\n",
    "The next step is to specify a forecasting metric. These are functions that return a number when input with prediction and actual series. They are different from `sklearn` metrics in that they accept series with indices rather than `np.array`s. Forecasting metrics can be invoked in two ways:\n",
    "\n",
    "* using the lean function interface, e.g., `mean_absolute_percentage_error` which is a python function `(y_true : pd.Series, y_pred : pd.Series) -> float`\n",
    "* using the composable class interface, e.g., `MeanAbsolutePercentageError`, which is a python class, callable with the same signature\n",
    "\n",
    "Casual users may opt to use the function interface. The class interface supports advanced use cases, such as parameter modification, custom metric composition, tuning over metric parameters (not covered in this tutorial)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.performance_metrics.forecasting import mean_absolute_percentage_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# option 1: using the lean function interface\n",
    "mean_absolute_percentage_error(y_test, y_pred)\n",
    "# note: the FIRST argument is the ground truth, the SECOND argument are the forecasts\n",
    "#       the order matters for most metrics in general"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To properly interpret numbers like this, it is useful to understand properties of the metric in question (e.g., lower is better), and to compare against suitable baselines and contender algorithms (see step 5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.performance_metrics.forecasting import MeanAbsolutePercentageError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# option 2: using the composable class interface\n",
    "mape = MeanAbsolutePercentageError(symmetric=False)\n",
    "# the class interface allows to easily construct variants of the MAPE\n",
    "#  e.g., the non-symmetric verion\n",
    "# it also allows for inspection of metric properties\n",
    "#  e.g., are higher values better (answer: no)?\n",
    "mape.greater_is_better"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluation works exactly like in option 2, but with the instantiated object\n",
    "mape(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: some metrics, such as `mean_absolute_scaled_error`, also require the training set for evaluation. In this case, the training set should be passed as a `y_train` argument. Refer to the API reference on individual metrics.\n",
    "\n",
    "NOTE: the workflow is the same for forecasters that make use of exogeneous data - no `X` is passed to the metrics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### step 5 - testing performance against benchmarks\n",
    "\n",
    "In general, forecast performances should be quantitatively tested against benchmark performances.\n",
    "\n",
    "Currently (`sktime` v0.6x), this is a roadmap development item. Contributions are very welcome."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.3.1 the basic batch forecast evaluation workflow in a nutshell - function metric interface<a class=\"anchor\" id=\"section_1_3_1\"></a>\n",
    "\n",
    "For convenience, we present the basic batch forecast evaluation workflow in one cell.\n",
    "This cell is using the lean function metric interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.base import ForecastingHorizon\n",
    "from sktime.forecasting.model_selection import temporal_train_test_split\n",
    "from sktime.forecasting.naive import NaiveForecaster\n",
    "from sktime.performance_metrics.forecasting import mean_absolute_percentage_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 1: splitting historical data\n",
    "y = load_airline()\n",
    "y_train, y_test = temporal_train_test_split(y, test_size=36)\n",
    "\n",
    "# step 2: running the basic forecasting workflow\n",
    "fh = ForecastingHorizon(y_test.index, is_relative=False)\n",
    "forecaster = NaiveForecaster(strategy=\"last\", sp=12)\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "\n",
    "# step 3: specifying the evaluation metric and\n",
    "# step 4: computing the forecast performance\n",
    "mean_absolute_percentage_error(y_test, y_pred)\n",
    "\n",
    "# step 5: testing forecast performance against baseline\n",
    "# under development"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.3.2 the basic batch forecast evaluation workflow in a nutshell - metric class interface<a class=\"anchor\" id=\"section_1_3_2\"></a>\n",
    "\n",
    "For convenience, we present the basic batch forecast evaluation workflow in one cell.\n",
    "This cell is using the advanced class specification interface for metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.base import ForecastingHorizon\n",
    "from sktime.forecasting.model_selection import temporal_train_test_split\n",
    "from sktime.forecasting.naive import NaiveForecaster\n",
    "from sktime.performance_metrics.forecasting import MeanAbsolutePercentageError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 1: splitting historical data\n",
    "y = load_airline()\n",
    "y_train, y_test = temporal_train_test_split(y, test_size=36)\n",
    "\n",
    "# step 2: running the basic forecasting workflow\n",
    "fh = ForecastingHorizon(y_test.index, is_relative=False)\n",
    "forecaster = NaiveForecaster(strategy=\"last\", sp=12)\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "\n",
    "# step 3: specifying the evaluation metric\n",
    "mape = MeanAbsolutePercentageError(symmetric=False)\n",
    "# if function interface is used, just use the function directly in step 4\n",
    "\n",
    "# step 4: computing the forecast performance\n",
    "mape(y_test, y_pred)\n",
    "\n",
    "# step 5: testing forecast performance against baseline\n",
    "# under development"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 advanced deployment workflow: rolling updates & forecasts<a class=\"anchor\" id=\"section_1_4\"></a>\n",
    "\n",
    "A common use case requires the forecaster to regularly update with new data and make forecasts on a rolling basis. This is especially useful if the same kind of forecast has to be made at regular time points, e.g., daily or weekly. `sktime` forecasters support this type of deployment workflow via the `update` and `update_predict` methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4.1 updating a forecaster with the `update` method<a class=\"anchor\" id=\"section_1_4_1\"></a>\n",
    "\n",
    "The `update` method can be called when a forecaster is already fitted, to ingest new data and make updated forecasts - this is referred to as an \"update step\".\n",
    "\n",
    "After the update, the forecaster's internal \"now\" state (the `cutoff`) is set to the latest time stamp seen in the update batch (assumed to be later than previously seen data).\n",
    "\n",
    "The general pattern is as follows:\n",
    "\n",
    "1. specify a forecasting strategy\n",
    "2. specify a relative forecasting horizon\n",
    "3. fit the forecaster to an initial batch of data using `fit`\n",
    "4. make forecasts for the relative forecasting horizon, using `predict`\n",
    "5. obtain new data; use `update` to ingest new data\n",
    "6. make forecasts using `predict` for the updated data\n",
    "7. repeat 5 and 6 as often as required\n",
    "\n",
    "**Example**: suppose that, in the airline example, we want to make forecasts a year ahead, but every month, starting December 1957. The first few months, forecasts would be made as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.ets import AutoETS\n",
    "from sktime.utils.plotting import plot_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we prepare the full data set for convenience\n",
    "# note that in the scenario we will \"know\" only part of this at certain time points\n",
    "y = load_airline()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# December 1957\n",
    "\n",
    "# this is the data known in December 1975\n",
    "y_1957Dec = y[:-36]\n",
    "\n",
    "# step 1: specifying the forecasting strategy\n",
    "forecaster = AutoETS(auto=True, sp=12, n_jobs=-1)\n",
    "\n",
    "# step 2: specifying the forecasting horizon: one year ahead, all months\n",
    "fh = np.arange(1, 13)\n",
    "\n",
    "# step 3: this is the first time we use the model, so we fit it\n",
    "forecaster.fit(y_1957Dec)\n",
    "\n",
    "# step 4: obtaining the first batch of forecasts for Jan 1958 - Dec 1958\n",
    "y_pred_1957Dec = forecaster.predict(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting predictions and past data\n",
    "plot_series(y_1957Dec, y_pred_1957Dec, labels=[\"y_1957Dec\", \"y_pred_1957Dec\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# January 1958\n",
    "\n",
    "# new data is observed:\n",
    "y_1958Jan = y[[-36]]\n",
    "\n",
    "# step 5: we update the forecaster with the new data\n",
    "forecaster.update(y_1958Jan)\n",
    "\n",
    "# step 6: making forecasts with the updated data\n",
    "y_pred_1958Jan = forecaster.predict(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note that the fh is relative, so forecasts are automatically for 1 month later\n",
    "#  i.e., from Feb 1958 to Jan 1959\n",
    "y_pred_1958Jan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting predictions and past data\n",
    "plot_series(\n",
    "    y[:-35],\n",
    "    y_pred_1957Dec,\n",
    "    y_pred_1958Jan,\n",
    "    labels=[\"y_1957Dec\", \"y_pred_1957Dec\", \"y_pred_1958Jan\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# February 1958\n",
    "\n",
    "# new data is observed:\n",
    "y_1958Feb = y[[-35]]\n",
    "\n",
    "# step 5: we update the forecaster with the new data\n",
    "forecaster.update(y_1958Feb)\n",
    "\n",
    "# step 6: making forecasts with the updated data\n",
    "y_pred_1958Feb = forecaster.predict(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting predictions and past data\n",
    "plot_series(\n",
    "    y[:-35],\n",
    "    y_pred_1957Dec,\n",
    "    y_pred_1958Jan,\n",
    "    y_pred_1958Feb,\n",
    "    labels=[\"y_1957Dec\", \"y_pred_1957Dec\", \"y_pred_1958Jan\", \"y_pred_1958Feb\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and so on.\n",
    "\n",
    "A shorthand for running first `update` and then `predict` is `update_predict_single` - for some algorithms, this may be more efficient than the separate calls to `update` and `predict`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# March 1958\n",
    "\n",
    "# new data is observed:\n",
    "y_1958Mar = y[[-34]]\n",
    "\n",
    "# step 5&6: update/predict in one step\n",
    "forecaster.update_predict_single(y_1958Mar, fh=fh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4.2 moving the \"now\" state without updating the model<a class=\"anchor\" id=\"section_1_4_2\"></a>\n",
    "\n",
    "In the rolling deployment mode, may be useful to move the estimator's \"now\" state (the `cutoff`) to later, for example if no new data was observed, but time has progressed; or, if computations take too long, and forecasts have to be queried.\n",
    "\n",
    "The `update` interface provides an option for this, via the `update_params` argument of `update` and other update funtions.\n",
    "\n",
    "If `update_params` is set to `False`, no model update computations are performed; only data is stored, and the internal \"now\" state (the `cutoff`) is set to the most recent date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# April 1958\n",
    "\n",
    "# new data is observed:\n",
    "y_1958Apr = y[[-33]]\n",
    "\n",
    "# step 5: perform an update without re-computing the model parameters\n",
    "forecaster.update(y_1958Apr, update_params=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4.3 walk-forward predictions on a batch of data<a class=\"anchor\" id=\"section_1_4_3\"></a>\n",
    "\n",
    "`sktime` can also simulate the update/predict deployment mode with a full batch of data.\n",
    "\n",
    "This is not useful in deployment, as it requires all data to be available in advance; however, it is useful in playback, such as for simulations or model evaluation.\n",
    "\n",
    "The update/predict playback mode can be called using `update_predict` and a re-sampling constructor which encodes the precise walk-forward scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from sktime.datasets import load_airline\n",
    "# from sktime.forecasting.ets import AutoETS\n",
    "# from sktime.forecasting.model_selection import ExpandingWindowSplitter\n",
    "# from sktime.utils.plotting import plot_series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: commented out - this part of the interface is currently undergoing a re-work. Contributions and PR are appreciated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for playback, the full data needs to be loaded in advance\n",
    "# y = load_airline()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 1: specifying the forecasting strategy\n",
    "# forecaster = AutoETS(auto=True, sp=12, n_jobs=-1)\n",
    "\n",
    "# step 2: specifying the forecasting horizon\n",
    "# fh - np.arange(1, 13)\n",
    "\n",
    "# step 3: specifying the cross-validation scheme\n",
    "# cv = ExpandingWindowSplitter()\n",
    "\n",
    "# step 4: fitting the forecaster - fh should be passed here\n",
    "# forecaster.fit(y[:-36], fh=fh)\n",
    "\n",
    "# step 5: rollback\n",
    "# y_preds = forecaster.update_predict(y, cv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 advanced evaluation worfklow: rolling re-sampling and aggregate errors, rolling back-testing<a class=\"anchor\" id=\"section_1_5\"></a>\n",
    "\n",
    "To evaluate forecasters with respect to their performance in rolling forecasting, the forecaster needs to be tested in a set-up mimicking rolling forecasting, usually on past data. Note that the batch back-testing as in Section 1.3 would not be an appropriate evaluation set-up for rolling deployment, as that tests only a single forecast batch. \n",
    "\n",
    "The advanced evaluation workflow can be carried out using the `evaluate` benchmarking function.\n",
    "`evalute` takes as arguments:\n",
    "- a `forecaster` to be evaluated\n",
    "- a `scikit-learn` re-sampling strategy for temporal splitting (`cv` below), e.g., `ExpandingWindowSplitter` or `SlidingWindowSplitter`\n",
    "- a `strategy` (string): whether the forecaster should be always be refitted or just fitted once and then updated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.arima import AutoARIMA\n",
    "from sktime.forecasting.model_evaluation import evaluate\n",
    "from sktime.forecasting.model_selection import ExpandingWindowSplitter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = AutoARIMA(sp=12, suppress_warnings=True)\n",
    "\n",
    "cv = ExpandingWindowSplitter(\n",
    "    step_length=12, fh=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], initial_window=72\n",
    ")\n",
    "\n",
    "df = evaluate(forecaster=forecaster, y=y, cv=cv, strategy=\"refit\", return_data=True)\n",
    "\n",
    "df.iloc[:, :5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualization of a forecaster evaluation\n",
    "fig, ax = plot_series(\n",
    "    y,\n",
    "    df[\"y_pred\"].iloc[0],\n",
    "    df[\"y_pred\"].iloc[1],\n",
    "    df[\"y_pred\"].iloc[2],\n",
    "    df[\"y_pred\"].iloc[3],\n",
    "    df[\"y_pred\"].iloc[4],\n",
    "    df[\"y_pred\"].iloc[5],\n",
    "    markers=[\"o\", \"\", \"\", \"\", \"\", \"\", \"\"],\n",
    "    labels=[\"y_true\"] + [\"y_pred (Backtest \" + str(x) + \")\" for x in range(6)],\n",
    ")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "todo: performance metrics, averages, and testing - contributions to `sktime` and the tutorial are welcome."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Forecasters in `sktime` - main families<a class=\"anchor\" id=\"chapter2\"></a>\n",
    "\n",
    "`sktime` supports a number of commonly used forecasters, many of them interfaced from state-of-art forecasting packages. All forecasters are available under the unified `sktime` interface.\n",
    "\n",
    "The main classes that are currently stably supported are:\n",
    "\n",
    "* `ExponentialSmoothing`, `ThetaForecaster`, and `autoETS` from `statsmodels`\n",
    "* `ARIMA` and `autoARIMA` from `pmdarima`\n",
    "* `BATS` and `TBATS` from `tbats`\n",
    "* `PolynomialTrend` for forecasting polynomial trends\n",
    "* `Prophet` which interfaces Facebook `prophet`\n",
    "\n",
    "For illustration, all estimators below will be presented on the basic forecasting workflow - though they also support the advanced forecasting and evaluation workflows under the unified `sktime` interface (see Section 1).\n",
    "\n",
    "For use in the other workflows, simply replace the \"forecaster specification block\" (\"`forecaster=`\") by the forecaster specification block in the examples presented below.\n",
    "\n",
    "Generally, all forecasters available in `sktime` can be listed with the `all_estimators` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.registry import all_estimators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all_estimators returns list of pairs - data frame conversion for pretty printing\n",
    "all_estimators(\"forecaster\", as_dataframe=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All forecasters follow the same interface, and can be used in the workflows presented in Section 1.\n",
    "\n",
    "We proceed by showcasing some commonnly used classes of forecasters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports necessary for this chapter\n",
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.base import ForecastingHorizon\n",
    "from sktime.forecasting.model_selection import temporal_train_test_split\n",
    "from sktime.performance_metrics.forecasting import mean_absolute_percentage_error\n",
    "from sktime.utils.plotting import plot_series\n",
    "\n",
    "# data loading for illustration (see section 1 for explanation)\n",
    "y = load_airline()\n",
    "y_train, y_test = temporal_train_test_split(y, test_size=36)\n",
    "fh = ForecastingHorizon(y_test.index, is_relative=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 exponential smoothing, theta forecaster, autoETS from `statsmodels`<a class=\"anchor\" id=\"section_2_1\"></a>\n",
    "\n",
    "`sktime` interfaces a number of statistical forecasting algorithms from `statsmodels`: exponential smoothing, theta, and aut-ETS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, to use exponential smoothing with an additive trend component and multiplicative seasonality on the airline data set, we can write the following.\n",
    "\n",
    "Note that since this is monthly data, a good choic for seasonal periodicity (sp) is 12 (= hypothesized periodicity of a year)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.exp_smoothing import ExponentialSmoothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = ExponentialSmoothing(trend=\"add\", seasonal=\"additive\", sp=12)\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exponential smoothing of state space model can also be automated similar\n",
    " to the [ets](https://www.rdocumentation.org/packages/forecast/versions/8.13/topics/ets) function in R. This is implemented in the `AutoETS` forecaster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.ets import AutoETS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = AutoETS(auto=True, sp=12, n_jobs=-1)\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# todo: explain Theta; explain how to get theta-lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 ARIMA and autoARIMA<a class=\"anchor\" id=\"section_2_2\"></a>\n",
    "\n",
    "`sktime` interfaces `pmdarima` for its ARIMA class models.\n",
    "For a classical ARIMA model with set parameters, use the `ARIMA` forecaster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.arima import ARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = ARIMA(\n",
    "    order=(1, 1, 0), seasonal_order=(0, 1, 0, 12), suppress_warnings=True\n",
    ")\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`AutoARIMA` is an automatically tuned `ARIMA` variant that obtains the optimal pdq parameters automatically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.arima import AutoARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = AutoARIMA(sp=12, suppress_warnings=True)\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to obtain the fitted parameters, run\n",
    "forecaster.get_fitted_params()\n",
    "# should these not include pdq?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 BATS and TBATS<a class=\"anchor\" id=\"section_2_3\"></a>\n",
    "\n",
    "`sktime` interfaces BATS and TBATS from the [`tbats`](https://github.com/intive-DataScience/tbats) package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.bats import BATS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = BATS(sp=12, use_trend=True, use_box_cox=False)\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.tbats import TBATS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = TBATS(sp=12, use_trend=True, use_box_cox=False)\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Facebook prophet<a class=\"anchor\" id=\"section_2_4\"></a>\n",
    "\n",
    "`sktime` provides an interface to [`fbprophet`](https://github.com/facebook/prophet) by Facebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.fbprophet import Prophet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The current interface does not support period indices, only pd.DatetimeIndex.\n",
    "\n",
    "Consider improving this by contributing the `sktime`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert index to pd.DatetimeIndex\n",
    "z = y.copy()\n",
    "z = z.to_timestamp(freq=\"M\")\n",
    "z_train, z_test = temporal_train_test_split(z, test_size=36)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = Prophet(\n",
    "    seasonality_mode=\"multiplicative\",\n",
    "    n_changepoints=int(len(y_train) / 12),\n",
    "    add_country_holidays={\"country_name\": \"Germany\"},\n",
    "    yearly_seasonality=True,\n",
    "    weekly_seasonality=False,\n",
    "    daily_seasonality=False,\n",
    ")\n",
    "\n",
    "forecaster.fit(z_train)\n",
    "y_pred = forecaster.predict(fh.to_relative(cutoff=y_train.index[-1]))\n",
    "y_pred.index = y_test.index\n",
    "\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 State Space Model (Structural Time Series)<a class=\"anchor\" id=\"section_2_5\"></a>\n",
    "\n",
    "We can also use the [`UobservedComponents`](https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.structural.UnobservedComponents.html) class from [`statsmodels`](https://www.statsmodels.org/stable/index.html) to generate predictions using a state space model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.structural import UnobservedComponents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can model seasonality using Fourier modes as in the Prophet model.\n",
    "forecaster = UnobservedComponents(\n",
    "    level=\"local linear trend\", freq_seasonal=[{\"period\": 12, \"harmonics\": 10}]\n",
    ")\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Advanced composition patterns - pipelines, reduction, autoML, and more<a class=\"anchor\" id=\"chapter3\"></a>\n",
    "\n",
    "`sktime` supports a number of advanced composition patterns to create forecasters out of simpler components:\n",
    "\n",
    "* reduction - building a forecaster from estimators of \"simpler\" scientific types, like `scikit-learn` regressors. A common example is feature/label tabulation by rolling window, aka the \"direct reduction strategy\".\n",
    "* tuning - determining values for hyper-parameters of a forecaster in a data-driven manner. A common example is grid search on temporally rolling re-sampling of train/test splits.\n",
    "* pipelining - concatenating transformers with a forecaster to obtain one forecaster. A common example is detrending and deseasonalizing then forecasting, an instance of this is the common \"STL forecaster\".\n",
    "* autoML, also known as automated model selection - using automated tuning strategies to select not only hyper-parameters but entire forecasting strategies. A common example is on-line multiplexer tuning.\n",
    "\n",
    "For illustration, all estimators below will be presented on the basic forecasting workflow - though they also support the advanced forecasting and evaluation workflows under the unified `sktime` interface (see Section 1).\n",
    "\n",
    "For use in the other workflows, simply replace the \"forecaster specification block\" (\"`forecaster=`\") by the forecaster specification block in the examples presented below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports necessary for this chapter\n",
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.base import ForecastingHorizon\n",
    "from sktime.forecasting.model_selection import temporal_train_test_split\n",
    "from sktime.performance_metrics.forecasting import mean_absolute_percentage_error\n",
    "from sktime.utils.plotting import plot_series\n",
    "\n",
    "# data loading for illustration (see section 1 for explanation)\n",
    "y = load_airline()\n",
    "y_train, y_test = temporal_train_test_split(y, test_size=36)\n",
    "fh = ForecastingHorizon(y_test.index, is_relative=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Reduction: from forecasting to regression<a class=\"anchor\" id=\"section_3_1\"></a>\n",
    "\n",
    "`sktime` provides a meta-estimator that allows the use of any `scikit-learn` estimator for forecasting.\n",
    "\n",
    "* **modular** and **compatible with scikit-learn**, so that we can easily apply any scikit-learn regressor to solve our forecasting problem,\n",
    "* **parametric** and **tuneable**, allowing us to tune hyper-parameters such as the window length or strategy to generate forecasts\n",
    "* **adaptive**, in the sense that it adapts the scikit-learn's estimator interface to that of a forecaster, making sure that we can tune and properly evaluate our model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**: we will define a tabulation reduction strategy to convert a k-nearest neighbors regressor (`sklearn` `KNeighborsRegressor`) into a forecaster. The composite algorithm is an object compliant with the `sktime` forecaster interface (picture: big robot), and contains the regressor as a parameter accessible component (picture: little robot). In `fit`, the composite algorithm uses a sliding window strategy to tabulate the data, and fit the regressor to the tabulated data (picture: left half). In `predict`, the composite algorithm presents the regressor with the last observed window to obtain predictions (picture: right half).\n",
    "\n",
    "<img src=\"img/forecasting-to-regression-reduction.png\" width=\"500\"/>\n",
    "\n",
    "Below, the composite is constructed using the shorthand function `make_reduction` which produces a `sktime` estimator of forecaster scitype. It is called with a constructed `scikit-learn` regressor, `regressor`, and additional parameter which can be later tuned as hyper-parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "\n",
    "from sktime.forecasting.compose import make_reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regressor = KNeighborsRegressor(n_neighbors=1)\n",
    "forecaster = make_reduction(regressor, window_length=15, strategy=\"recursive\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:07.439875Z",
     "iopub.status.busy": "2021-04-10T16:07:07.439381Z",
     "iopub.status.idle": "2021-04-10T16:07:07.650204Z",
     "shell.execute_reply": "2021-04-10T16:07:07.650693Z"
    }
   },
   "outputs": [],
   "source": [
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example we use the \"recursive\" reduction strategy. Other implemented strategies are: \n",
    "* \"direct\", \n",
    "* \"dirrec\", \n",
    "* \"multioutput\". "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters can be inspected using `scikit-learn` compatible `get_params` functionality (and set using `set_params`). This provides tunable and nested access to parameters of the `KNeighborsRegressor` (as `estimator_etc`), and the `window_length` of the reduction strategy. Note that the `strategy` is not accessible, as underneath the utility function this is mapped on separate algorithm classes. For tuning over algorithms, see the \"autoML\" section below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster.get_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Pipelining, detrending and deseasonalization<a class=\"anchor\" id=\"section_3_2\"></a>\n",
    "\n",
    "A common composition motif is pipelining: for example, first deseasonalizing or detrending the data, then forecasting the detrended/deseasonalized series. When forecasting, one needs to add the trend and seasonal component back to the data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.1 The basic forecasting pipeline<a class=\"anchor\" id=\"section_3_2_1\"></a>\n",
    "\n",
    "`sktime` provides a generic pipeline object for this kind of composite modelling, the `TransforemdTargetForecaster`. It chains an arbitrary number of transformations with a forecaster. The transformations should be instances of estimators with series-to-series-transformer scitype. An example of the syntax is below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.arima import ARIMA\n",
    "from sktime.forecasting.compose import TransformedTargetForecaster\n",
    "from sktime.transformations.series.detrend import Deseasonalizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = TransformedTargetForecaster(\n",
    "    [\n",
    "        (\"deseasonalize\", Deseasonalizer(model=\"multiplicative\", sp=12)),\n",
    "        (\"forecast\", ARIMA()),\n",
    "    ]\n",
    ")\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `TransformedTargetForecaster` is constructed with a list of steps, each a pair of name and estimator. The last estimator should be of forecaster scitype, the other estimators should be series-to-series transformers which possess both a `transform` and `inverse_transform` method. The resulting estimator is of forecaster scitype and has all interface defining methods. In `fit`, all transformers apply `fit_transforms` to the data, then the forecaster's `fit`; in `predict`, first the forecaster's `predict` is applied, then the transformers' `inverse_transform` in reverse order."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.2 The `Detrender` as pipeline component<a class=\"anchor\" id=\"section_3_2_2\"></a>\n",
    "\n",
    "For detrending, we can use the `Detrender`. This is an estimator of series-to-transformer scitype that wraps an arbitrary forecaster. For example, for linear detrending, we can use `PolynomialTrendForecaster` to fit a linear trend, and then subtract/add it using the `Detrender` transformer inside `TransformedTargetForecaster`.\n",
    "\n",
    "To understand better what happens, we first examine the detrender separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.trend import PolynomialTrendForecaster\n",
    "from sktime.transformations.series.detrend import Detrender"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# linear detrending\n",
    "forecaster = PolynomialTrendForecaster(degree=1)\n",
    "transformer = Detrender(forecaster=forecaster)\n",
    "yt = transformer.fit_transform(y_train)\n",
    "\n",
    "# internally, the Detrender uses the in-sample predictions\n",
    "# of the PolynomialTrendForecaster\n",
    "forecaster = PolynomialTrendForecaster(degree=1)\n",
    "fh_ins = -np.arange(len(y_train))  # in-sample forecasting horizon\n",
    "y_pred = forecaster.fit(y_train).predict(fh=fh_ins)\n",
    "\n",
    "plot_series(y_train, y_pred, yt, labels=[\"y_train\", \"fitted linear trend\", \"residuals\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the `Detrender` is of scitype series-to-series-transformer, it can be used in the `TransformedTargetForecaster` for detrending any forecaster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = TransformedTargetForecaster(\n",
    "    [\n",
    "        (\"deseasonalize\", Deseasonalizer(model=\"multiplicative\", sp=12)),\n",
    "        (\"detrend\", Detrender(forecaster=PolynomialTrendForecaster(degree=1))),\n",
    "        (\"forecast\", ARIMA()),\n",
    "    ]\n",
    ")\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.3 Complex pipeline composites and parameter inspection<a class=\"anchor\" id=\"section_3_2_3\"></a>\n",
    "\n",
    "`sktime` follows the `scikit-learn` philosophy of composability and nested parameter inspection. As long as an estimator has the right scitype, it can be used as part of any composition principle requiring that scitype. Above, we have already seen the example of a forecaster inside a `Detrender`, which is an estimator of scitype series-to-series-transformer, with one component of forecaster scitype. Similarly, in a `TransformedTargetForecaster`, we can use the reduction composite from Section 3.1 as the last forecaster element in the pipeline, which inside has an estimator of tabular regressor scitype, the `KNeighborsRegressor`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "\n",
    "from sktime.forecasting.compose import make_reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = TransformedTargetForecaster(\n",
    "    [\n",
    "        (\"deseasonalize\", Deseasonalizer(model=\"multiplicative\", sp=12)),\n",
    "        (\"detrend\", Detrender(forecaster=PolynomialTrendForecaster(degree=1))),\n",
    "        (\n",
    "            \"forecast\",\n",
    "            make_reduction(\n",
    "                KNeighborsRegressor(),\n",
    "                scitype=\"tabular-regressor\",\n",
    "                window_length=15,\n",
    "                strategy=\"recursive\",\n",
    "            ),\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with `scikit-learn` models, we can inspect and access parameters  of any component via `get_params` and `set_params`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster.get_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Parameter tuning<a class=\"anchor\" id=\"section_3_3\"></a>\n",
    "\n",
    "`sktime` provides parameter tuning strategies as compositors of forecaster scitype, similar to `scikit-learn`'s `GridSearchCV`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3.1 Basic tuning using `ForecastingGridSearchCV`<a class=\"anchor\" id=\"section_3_3_1\"></a>\n",
    "\n",
    "The compositor `ForecastingGridSearchCV` (and other tuners) are constructed with a forecaster to tune, a cross-validation constructor, a `scikit-learn` parameter grid, and parameters specific to the tuning strategy. Cross-validation constructors follow the `scikit-learn` interface for re-samplers, and can be slotted in exchangeably.\n",
    "\n",
    "As an example, we show tuning of the window length in the reduction compositor from Section 3.1, using temporal sliding window tuning:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "\n",
    "from sktime.forecasting.compose import make_reduction\n",
    "from sktime.forecasting.model_selection import (\n",
    "    ForecastingGridSearchCV,\n",
    "    SlidingWindowSplitter,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:36.557083Z",
     "iopub.status.busy": "2021-04-10T16:07:36.495720Z",
     "iopub.status.idle": "2021-04-10T16:07:37.096548Z",
     "shell.execute_reply": "2021-04-10T16:07:37.097150Z"
    }
   },
   "outputs": [],
   "source": [
    "regressor = KNeighborsRegressor()\n",
    "forecaster = make_reduction(regressor, window_length=15, strategy=\"recursive\")\n",
    "param_grid = {\"window_length\": [7, 12, 15]}\n",
    "\n",
    "# We fit the forecaster on an initial window which is 80% of the historical data\n",
    "# then use temporal sliding window cross-validation to find the optimal hyper=parameters\n",
    "cv = SlidingWindowSplitter(initial_window=int(len(y_train) * 0.8), window_length=20)\n",
    "gscv = ForecastingGridSearchCV(\n",
    "    forecaster, strategy=\"refit\", cv=cv, param_grid=param_grid\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with other composites, the resulting forecaster provides the unified interface of `sktime` forecasters - window splitting, tuning, etc requires no manual effort and is done behind the unified interface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:37.120099Z",
     "iopub.status.busy": "2021-04-10T16:07:37.100367Z",
     "iopub.status.idle": "2021-04-10T16:07:37.372802Z",
     "shell.execute_reply": "2021-04-10T16:07:37.372246Z"
    }
   },
   "outputs": [],
   "source": [
    "gscv.fit(y_train)\n",
    "y_pred = gscv.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tuned parameters can be accessed in the `best_params_` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:37.375928Z",
     "iopub.status.busy": "2021-04-10T16:07:37.375391Z",
     "iopub.status.idle": "2021-04-10T16:07:37.377362Z",
     "shell.execute_reply": "2021-04-10T16:07:37.377899Z"
    }
   },
   "outputs": [],
   "source": [
    "gscv.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An instance of the best forecaster, with hyper-parameters set, can be retrieved by accessing the `best_forecaster_` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.best_forecaster_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3.2 Tuning of complex composites<a class=\"anchor\" id=\"section_3_3_2\"></a>\n",
    "\n",
    "As in `scikit-learn`, parameters of nested components can be tuned by accessing their `get_params` key - by default this is `[estimatorname]__[parametername]` if `[estimatorname]` is the name of the component, and `[parametername]` the name of a parameter within the estimator `[estimatorname]`. \n",
    "\n",
    "For example, below we tune the `KNeighborsRegressor` component's `n_neighbors`, in addition to tuning `window_length`. The tuneable parameters can easily be queried using `forecaster.get_params()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "\n",
    "from sktime.forecasting.compose import make_reduction\n",
    "from sktime.forecasting.model_selection import (\n",
    "    ForecastingGridSearchCV,\n",
    "    SlidingWindowSplitter,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "param_grid = {\"window_length\": [7, 12, 15], \"estimator__n_neighbors\": np.arange(1, 10)}\n",
    "\n",
    "regressor = KNeighborsRegressor()\n",
    "forecaster = make_reduction(\n",
    "    regressor, scitype=\"tabular-regressor\", strategy=\"recursive\"\n",
    ")\n",
    "\n",
    "cv = SlidingWindowSplitter(initial_window=int(len(y_train) * 0.8), window_length=30)\n",
    "gscv = ForecastingGridSearchCV(forecaster, cv=cv, param_grid=param_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.fit(y_train)\n",
    "y_pred = gscv.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative to the above is tuning the regressor separately, using `scikit-learn`'s `GridSearchCV` and a separate parameter grid. As this does not use the \"overall\" performance metric to tune the inner regressor, performance of the composite forecaster may vary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "# tuning the 'n_estimator' hyperparameter of RandomForestRegressor from scikit-learn\n",
    "regressor_param_grid = {\"n_neighbors\": np.arange(1, 10)}\n",
    "forecaster_param_grid = {\"window_length\": [7, 12, 15]}\n",
    "\n",
    "# create a tunnable regressor with GridSearchCV\n",
    "regressor = GridSearchCV(KNeighborsRegressor(), param_grid=regressor_param_grid)\n",
    "forecaster = make_reduction(\n",
    "    regressor, scitype=\"tabular-regressor\", strategy=\"recursive\"\n",
    ")\n",
    "\n",
    "cv = SlidingWindowSplitter(initial_window=int(len(y_train) * 0.8), window_length=30)\n",
    "gscv = ForecastingGridSearchCV(forecaster, cv=cv, param_grid=forecaster_param_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.fit(y_train)\n",
    "y_pred = gscv.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: a smart implementation of this would use caching to save partial results from the inner tuning and reduce runtime substantially - currently `sktime` does not support this. Consider helping to improve `sktime`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3.3 Selecting the metric and retrieving scores<a class=\"anchor\" id=\"section_3_3_3\"></a>\n",
    "\n",
    "All tuning algorithms in `sktime` allow the user to set a score; for forecasting the default is mean absolute percentage error. The score can be set using the `score` argument, to any scorer function or class, as in Section 1.3.\n",
    "\n",
    "Re-sampling tuners retain performances on individual forecast re-sample folds, which can be retrieved from the `cv_results_` argument after the forecaster has been fit via a call to `fit`.\n",
    "\n",
    "In the above example, using the mean squared error instead of the mean absolute percentage error for tuning would be done by defining the forecaster as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.performance_metrics.forecasting import MeanSquaredError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mse = MeanSquaredError()\n",
    "\n",
    "param_grid = {\"window_length\": [7, 12, 15]}\n",
    "\n",
    "regressor = KNeighborsRegressor()\n",
    "cv = SlidingWindowSplitter(initial_window=int(len(y_train) * 0.8), window_length=30)\n",
    "\n",
    "gscv = ForecastingGridSearchCV(forecaster, cv=cv, param_grid=param_grid, scoring=mse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The performances on individual folds can be accessed as follows, after fitting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.fit(y_train)\n",
    "gscv.cv_results_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  3.4 autoML  aka automated model selection, ensembling and hedging<a class=\"anchor\" id=\"section_3_4\"></a>\n",
    "\n",
    "`sktime` provides a number of compositors for ensembling and automated model selection. In contrast to tuning, which uses data-driven strategies to find optimal hyper-parameters for a fixed forecaster, the strategies in this section combine or select on the level of estimators, using a collection of forecasters to combine or select from.\n",
    "\n",
    "The strategies discussed in this section are:\n",
    "* autoML aka automated model selection\n",
    "* simple ensembling\n",
    "* prediction weighted ensembles with weight updates, and hedging strategies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  3.4.1 autoML aka automatic model selection, using tuning plus multiplexer<a class=\"anchor\" id=\"section_3_4_1\"></a>\n",
    "\n",
    "The most flexible way to perform model selection over forecasters is by using the `MultiplexForecaster`, which exposes the choice of a forecaster from a list as a hyper-parameter that is tunable by generic hyper-parameter tuning strategies such as in Section 3.3.\n",
    "\n",
    "In isolation, `MultiplexForecaster` is constructed with a named list `forecasters`, of forecasters. It has a single hyper-parameter, `selected_forecaster`, which can be set to the name of any forecaster in `forecasters`, and behaves exactly like the forecaster keyed in `forecasters` by `selected_forecaster`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.compose import MultiplexForecaster\n",
    "from sktime.forecasting.exp_smoothing import ExponentialSmoothing\n",
    "from sktime.forecasting.naive import NaiveForecaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:14:33.332255Z",
     "iopub.status.busy": "2021-04-10T16:14:33.310571Z",
     "iopub.status.idle": "2021-04-10T16:14:34.759567Z",
     "shell.execute_reply": "2021-04-10T16:14:34.760081Z"
    }
   },
   "outputs": [],
   "source": [
    "forecaster = MultiplexForecaster(\n",
    "    forecasters=[\n",
    "        (\"naive\", NaiveForecaster(strategy=\"last\")),\n",
    "        (\"ets\", ExponentialSmoothing(trend=\"add\", sp=12)),\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster.set_params(**{\"selected_forecaster\": \"naive\"})\n",
    "# now forecaster behaves like NaiveForecaster(strategy=\"last\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster.set_params(**{\"selected_forecaster\": \"ets\"})\n",
    "# now forecaster behaves like ExponentialSmoothing(trend=\"add\", sp=12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `MultiplexForecaster` is not too useful in isolation, but allows for flexible autoML when combined with a tuning wrapper. The below defines a forecaster that selects one of `NaiveForecaster` and `ExponentialSmoothing` by sliding window tuning as in Section 3.3.\n",
    "\n",
    "Combined with rolling use of the forecaster via the `update` functionality (see Section 1.4), the tuned multiplexer can switch back and forth between `NaiveForecaster` and `ExponentialSmoothing`, depending on performance, as time progresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.model_selection import (\n",
    "    ForecastingGridSearchCV,\n",
    "    SlidingWindowSplitter,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecaster = MultiplexForecaster(\n",
    "    forecasters=[\n",
    "        (\"naive\", NaiveForecaster(strategy=\"last\")),\n",
    "        (\"ets\", ExponentialSmoothing(trend=\"add\", sp=12)),\n",
    "    ]\n",
    ")\n",
    "cv = SlidingWindowSplitter(initial_window=int(len(y_train) * 0.5), window_length=30)\n",
    "forecaster_param_grid = {\"selected_forecaster\": [\"ets\", \"naive\"]}\n",
    "gscv = ForecastingGridSearchCV(forecaster, cv=cv, param_grid=forecaster_param_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.fit(y_train)\n",
    "y_pred = gscv.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with any tuned forecaster, best parameters and an instance of the tuned forecaster can be retrieved using `best_params_` and `best_forecaster_`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:14:35.007827Z",
     "iopub.status.busy": "2021-04-10T16:14:35.007334Z",
     "iopub.status.idle": "2021-04-10T16:14:35.009284Z",
     "shell.execute_reply": "2021-04-10T16:14:35.009782Z"
    }
   },
   "outputs": [],
   "source": [
    "gscv.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.best_forecaster_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4.2 autoML: selecting transformer combinations via `OptimalPassthrough`<a class=\"anchor\" id=\"section_3_4_2\"></a>\n",
    "\n",
    "`sktime` also provides capabilities for automated selection of pipeline components *inside* a pipeline, i.e., pipeline structure. This is achieved with the `OptionalPassthrough` transformer.\n",
    "\n",
    "The `OptionalPassthrough` transformer allows to tune whether a transformer inside a pipeline is applied to the data or not. For example, if we want to tune whether `sklearn.StandardScaler` is bringing an advantage to the forecast or not, we wrap it in `OptionalPassthrough`. Internally, `OptionalPassthrough` has a hyperparameter `passthrough: bool` that is tuneable; when `False` the composite behaves like the wrapped transformer, when `True`, it ignores the transformer within.\n",
    "\n",
    "To make effective use of `OptionalPasstrhough`, define a suitable parameter set using the `__` (double underscore) notation familiar from `scikit-learn`. This allows to access and tune attributes of nested objects like TabularToSeriesAdaptor(StandardScaler()). We can use `__` multiple times if we have more than two levels of nesting.\n",
    "\n",
    "In the following example, we take a deseasonalize/scale pipeline and tune over the four possible combinations of deseasonalizer and scaler being included in the pipeline yes/no (2 times 2 = 4); as well as over the forecaster's and the scaler's parameters.\n",
    "\n",
    "Note: this could be arbitrarily combined with `MultiplexForecaster`, as in Section 3.4.1, to select over pipeline architecture as well as over pipeline structure.\n",
    "\n",
    "Note: `scikit-learn` and `sktime` do not support conditional parameter sets at current (unlike, e.g., the `mlr3` package). This means that the grid search will optimize over the `scaler`'s parameters even when it is skipped. Designing/implementing this capability would be an interesting area for contributions or research."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "from sktime.datasets import load_airline\n",
    "from sktime.forecasting.compose import TransformedTargetForecaster\n",
    "from sktime.forecasting.model_selection import (\n",
    "    ForecastingGridSearchCV,\n",
    "    SlidingWindowSplitter,\n",
    ")\n",
    "from sktime.forecasting.naive import NaiveForecaster\n",
    "from sktime.transformations.series.adapt import TabularToSeriesAdaptor\n",
    "from sktime.transformations.series.compose import OptionalPassthrough\n",
    "from sktime.transformations.series.detrend import Deseasonalizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create pipeline\n",
    "pipe = TransformedTargetForecaster(\n",
    "    steps=[\n",
    "        (\"deseasonalizer\", OptionalPassthrough(Deseasonalizer())),\n",
    "        (\"scaler\", OptionalPassthrough(TabularToSeriesAdaptor(StandardScaler()))),\n",
    "        (\"forecaster\", NaiveForecaster()),\n",
    "    ]\n",
    ")\n",
    "\n",
    "# putting it all together in a grid search\n",
    "cv = SlidingWindowSplitter(\n",
    "    initial_window=60, window_length=24, start_with_window=True, step_length=24\n",
    ")\n",
    "param_grid = {\n",
    "    \"deseasonalizer__passthrough\": [True, False],\n",
    "    \"scaler__transformer__transformer__with_mean\": [True, False],\n",
    "    \"scaler__passthrough\": [True, False],\n",
    "    \"forecaster__strategy\": [\"drift\", \"mean\", \"last\"],\n",
    "}\n",
    "gscv = ForecastingGridSearchCV(forecaster=pipe, param_grid=param_grid, cv=cv, n_jobs=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gscv.fit(y_train)\n",
    "y_pred = gscv.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  3.4.3 simple ensembling strategies<a class=\"anchor\" id=\"section_3_4_3\"></a>\n",
    "\n",
    "TODO - contributions in this section are appreciated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sktime.forecasting.compose import EnsembleForecaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:07:36.322267Z",
     "iopub.status.busy": "2021-04-10T16:07:36.287424Z",
     "iopub.status.idle": "2021-04-10T16:07:36.486587Z",
     "shell.execute_reply": "2021-04-10T16:07:36.487094Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ses = ExponentialSmoothing(sp=12)\n",
    "holt = ExponentialSmoothing(trend=\"add\", damped_trend=False, sp=12)\n",
    "damped = ExponentialSmoothing(trend=\"add\", damped_trend=True, sp=12)\n",
    "\n",
    "forecaster = EnsembleForecaster(\n",
    "    [\n",
    "        (\"ses\", ses),\n",
    "        (\"holt\", holt),\n",
    "        (\"damped\", damped),\n",
    "    ]\n",
    ")\n",
    "forecaster.fit(y_train)\n",
    "y_pred = forecaster.predict(fh)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4.4 Prediction weighted ensembles and hedge ensembles<a class=\"anchor\" id=\"section_3_4_4\"></a>\n",
    "\n",
    "For model evaluation, we sometimes want to evaluate multiple forecasts, using temporal cross-validation with a sliding window over the test data. For this purpose, we can leverage the forecasters from the `online_forecasting` module which use a composite forecaster, `PredictionWeightedEnsemble`, to keep track of the loss accumulated by each forecaster and create a prediction weighted by the predictions of the most \"accurate\" forecasters.\n",
    "\n",
    "Note that the forecasting task is changed: we make 35 predictions since we need the first prediction to help update the weights, we do not predict 36 steps ahead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:14:39.872044Z",
     "iopub.status.busy": "2021-04-10T16:14:39.871304Z",
     "iopub.status.idle": "2021-04-10T16:14:39.874814Z",
     "shell.execute_reply": "2021-04-10T16:14:39.875324Z"
    }
   },
   "outputs": [],
   "source": [
    "from sktime.forecasting.all import mean_squared_error\n",
    "from sktime.forecasting.online_learning import (\n",
    "    NormalHedgeEnsemble,\n",
    "    OnlineEnsembleForecaster,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to initialize a `PredictionWeightedEnsembler` that will keep track of the loss accumulated by each forecaster and define which loss function we would like to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:14:39.878435Z",
     "iopub.status.busy": "2021-04-10T16:14:39.877913Z",
     "iopub.status.idle": "2021-04-10T16:14:39.879324Z",
     "shell.execute_reply": "2021-04-10T16:14:39.879979Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "hedge_expert = NormalHedgeEnsemble(n_estimators=3, loss_func=mean_squared_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We can then create the forecaster by defining the individual forecasters and specifying the `PredictionWeightedEnsembler` we are using. Then by fitting our forecasters and performing updates and prediction with the `update_predict` function, we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-10T16:14:39.923246Z",
     "iopub.status.busy": "2021-04-10T16:14:39.920259Z",
     "iopub.status.idle": "2021-04-10T16:14:40.917504Z",
     "shell.execute_reply": "2021-04-10T16:14:40.918060Z"
    }
   },
   "outputs": [],
   "source": [
    "forecaster = OnlineEnsembleForecaster(\n",
    "    [\n",
    "        (\"ses\", ses),\n",
    "        (\"holt\", holt),\n",
    "        (\"damped\", damped),\n",
    "    ],\n",
    "    ensemble_algorithm=hedge_expert,\n",
    ")\n",
    "\n",
    "forecaster.fit(y=y_train, fh=fh)\n",
    "y_pred = forecaster.update_predict_single(y_test)\n",
    "plot_series(y_train, y_test, y_pred, labels=[\"y_train\", \"y_test\", \"y_pred\"])\n",
    "mean_absolute_percentage_error(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Extension guide - implementing your own forecaster<a class=\"anchor\" id=\"chapter4\"></a>\n",
    "\n",
    "`sktime` is meant to be easily extensible, for direct contribution to `sktime` as well as for local/private extension with custom methods.\n",
    "\n",
    "To extend `sktime` with a new local or contributed forecaster, a good workflow to follow is:\n",
    "\n",
    "1. read through the [forecasting extension template](https://github.com/alan-turing-institute/sktime/blob/main/extension_templates/forecasting.py) - this is a `python` file with `todo` blocks that mark the places in which changes need to be added.\n",
    "2. optionally, if you are planning any major surgeries to the interface: look at the [base class architecture](https://github.com/alan-turing-institute/sktime/blob/main/sktime/forecasting/base/_base.py) - note that \"ordinary\" extension (e.g., new algorithm) should be easily doable without this.\n",
    "3. copy the forecasting extension template to a local folder in your own repository (local/private extension), or to a suitable location in your clone of the `sktime` or affiliated repository (if contributed extension), inside `sktime.forecasting`; rename the file and update the file docstring appropriately.\n",
    "4. address the \"todo\" parts. Usually, this means: changing the name of the class, setting the tag values, specifying hyper-parameters, filling in `__init__`, `_fit`, `_predict`, and optional methods such as `_update` (for details see the extension template). You can add private methods as long as they do not override the default public interface. For more details, see the extension template.\n",
    "5. to test your estimator manually: import your estimator and run it in the worfklows in Section 1; then use it in the compositors in Section 3.\n",
    "6. to test your estimator automatically: call `sktime.tests.test_all_estimators.test_estimator` on your estimator - note that the function takes the class, not an object instance. Before the call, you need to register the new estimator in `sktime.tests._config`, as an import, and by adding default parameter settings to the `ESTIMATOR_TEST_PARAMS` variable (the `dict` entry `key` is the class, and entry is a `scikit-learn` parameter set). `pytest` will also add the call to its automated tests in a working clone of the `sktime` repository.\n",
    "\n",
    "In case of direct contribution to `sktime` or one of its affiliated packages, additionally:\n",
    "* add yourself as an author to the code, and to the `CODEOWNERS` for the new estimator file(s).\n",
    "* create a pull request that contains only the new estimators (and their inheritance tree, if it's not just one class), as well as the automated tests as described above.\n",
    "* in the pull request, describe the estimator and optimally provide a publication or other technical reference for the strategy it implements.\n",
    "* before making the pull request, ensure that you have all necessary permissions to contribute the code to a permissive license (BSD-3) open source project."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Summary<a class=\"anchor\" id=\"chapter5\"></a>\n",
    "\n",
    "* `sktime` comes with several forecasting algorithms (or forecasters), all of which share a common interface. The interface is fully interoperable with the `scikit-learn` interface, and provides dedicated interface points for forecasting in batch and rolling mode.\n",
    "\n",
    "* `sktime` comes with rich composition functionality that allows to build complex pipelines easily, and connect easily with other parts of the open source ecosystem, such as `scikit-learn` and individual algorithm libraries.\n",
    "\n",
    "* `sktime` is easy to extend, and comes with user friendly tools to facilitate implementing and testing your own forecasters and composition principles.\n",
    "\n",
    "\n",
    "## Useful resources\n",
    "* For more details, take a look at [our paper on forecasting with sktime](https://arxiv.org/abs/2005.08067) in which we discuss the forecasting API in more detail and use it to replicate and extend the M4 study.\n",
    "* For a good introduction to forecasting, see [Hyndman, Rob J., and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2018](https://otexts.com/fpp2/).\n",
    "* For comparative benchmarking studies/forecasting competitions, see the [M4 competition](https://www.sciencedirect.com/science/article/pii/S0169207019301128) and the [M5 competition](https://www.kaggle.com/c/m5-forecasting-accuracy/overview)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "interpreter": {
   "hash": "fcc5fed35031463a248402718f3bbb1a61c709e60741a9777b2268658fc045fd"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
